Correlations of quasi-2D dipolar ultracold gas at finite temperatures 
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We study a quasi two dimensional dipolar gas at finite, but ultralow temperatures using the classical field 
approximation. The method, already used for a contact interacting gas, is extended here to samples with a 
weakly interacting long-range inter-atomic potential. We present statistical properties of the system for the 
current experiment with Chromium [1] and compare them with statistics for atoms with larger magnetic dipole 
moments. Significant enhancement of the third order correlation function, relevant for the particle losses, is 
found. 
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The successful cooling down of dipolar gases below the 
condensation temperature [2] has attracted many theorists and 
experimental groups. Dipolar forces may stabilize or destabi- 
lize the ultracold cloud depending on its polarization. A col- 
lapsing cloud forms so called Bose-Nova [3]. Dipolar forces 
also introduce inter-site effects between atoms loaded in opti- 
cal lattices [1]. From the theory side, a variety of new phe- 
nomena are expected, i.e. new quantum phases [4] or ro- 
tonization of the spectrum [5]. The interest of the community 
is still increasing, as the new atomic species. Erbium [6] and 
Dysprosium [7], has been condensed recently. Furthermore, 
a big effort has been made to condense polar molecules with 
electric dipole interactions orders of magnitude larger than the 
magnetic ones [8, 9]. 

Especially interesting is a quasi two dimensional system, 
where roton may be formed [10]. The latter, although investi- 
gated in many theoretical papers, has not been observed in 
dipolar gases yet. On the other hand it is known, that the 
two dimensional system at finite temperatures has many in- 
teresting features even in a purely contact interacting gas. The 
most intriguing is the Berezinsky-Kosterlitz-Thouless transi- 
tion, also studied experimentally [11]. 

Thus, it is natural to look at a two dimensional dipolar gas 
at temperature above the absolute zero. The method applied 
most often for such study is the Bogoliubov approximation 
[10]. Such treatment neglects higher order interactions be- 
tween modes and is limited to low temperatures only. More- 
over in the case of an untrapped system, it assumes a con- 
densation in the zero-momentum mode, whereas close to in- 
stability this is not necessarily the case. For system in a trap 
the appropriate shape of the condensate and, partially, interac- 
tions between excitations have been included in extensions of 
the Bogoliubov approximation [12-14], however, at the cost 
of more computationally expensive algorithms. Another qual- 
itative approach is based on the mean field method [15]. 

Among other very promising techniques are the classical 
field methods (called also c- field methods). In the case of a 
trapped system the numerical algorithm for a projected Gross- 
Pitaevskii equation has already been presented [16]. 

We present an even simpler example of the c-field meth- 
ods' family. It has been successfully used for deriving the full 
statistics in a thermal equilibrium of an ideal gas [17, 18] and 
a contact interacting gas for weakly repulsive [19, 20], as well 
as attractive forces [21]. Here, we apply this method to the 



dipolar gases, trapped in a "hybrid" geometry - a strong har- 
monic trap in the z direction and a box with periodic boundary 
condition in the perpendicular directions. With our method we 
gain qualitative insight into the statistical properties of the sys- 
tem, without any ad-hoc assumption for a shape of the mostly 
occupied state and including interactions in all orders. We 
use the canonical ensemble, thus the total number of atoms is 
conserved. 

The paper is organized as follow: In Sec. I we sketch our 
realization of the classical field approximation (CFA). The 
derivation of the dipolar interaction energy functionals, has 
been relegated to the Appendix. Sec. II focuses on the analy- 
sis of the, up to date, the best experimental realization of quasi 
2D dipolar gases, done with Chromium in [1], and compare it 
with Erbium and Dysprosium. The most important results are 
summarized in Conclusions III. 



I. CLASSICAL FIELD APPROXIMATION 

We study a weakly interacting dipolar Bose gas trapped in 
a very steep harmonic potential in the Z direction and square 
box with periodic boundary condition in the X and Y direc- 
tions. Thus, our Hamiltonian of the three dimensional Bose 
gas has a form: 
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+ i f d^rd^r'¥{r)^''{r')Vd^pir-r')il{r')^{r). 

The Hamiltonian is a sum of a single particle energy, a conven- 
tional contact interaction energy with the coupling constant g 
and a dipolar interaction energy. The dipolar potential is given 
by: 



Vdip (r) 



^dd r^ - 3 (nr) 



(2) 



where n is a unit vector in the direction of polarization and gdd 
parameterizes strength of dipolar interactions. The periodic 
boundary conditions in the X, Y directions are responsible 



for the quantization of momenta: 



i^i — J \}x-!'^y\i ^ — [f^it ^z\ 



(3) 



where L is the size of the box and ix, iy are any integer num- 
bers. Here we use the symbol i to denote the two dimensional 
index. The two dimensional vector fcj is a projection of a total 
momentum k onto Xy-plane, whereas k^ is the Z component 
of k. We restrict our considerations to such strong frequencies 
uiz, small interaction energies and low temperatures, that the 
shape of the atomic cloud does not change in the Z direction. 
In other words we assume, that a convenient base is provided 
by the product of the ground state Gaussian function in the 
strongly confined Z direction 
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with Iz being the Gaussian width (dispersion), and plane 
waves in other directions. Thus the suitable expansion of the 
field operator is given by: 
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where aj annihilates an atom in the i-th momentum state, and 
r = (p,z). 

The classical field approximation consists of replacing the 
creation and annihilation operators by classical complex am- 
plitudes: 
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It has been shown in [17], that probabilistic properties of the 
condensate for two dimensional ideal Bose gas in a box are for 
large N perfectly reproduced by the classical fields approxi- 
mation provided the number of degrees of freedom is kept 
finite with the last retained state chosen as: 



rfc^„^/2m = O.eSfcsT 
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where T is the absolute temperature and ks is the Boltzmann 
constant. However for small systems, the formula (7) needs a 
significant correction. In general, the optimal cut-off momen- 
tum kjnax is obtained from fitting the distribution of the num- 
ber of condensed atoms from the classical field to the exact 
solution in the ideal gas case. The latter is obtained according 
to [22], the former can be also exactly calculated as described 
in the Appendix of [19]. For instance, in the system described 
in the next section, the appropriate cut-off is given by: 



h'^kl^^j2m = 0.99fcsT, 



(8) 



a result which is obtained by a direct calculation. This ver- 
sion of classical fields approach was successfully used for the 
weakly contact interacting gas [18, 19, 21]. It also provided a 
quantitative description of recent measurements (supplemen- 
tal material in [20]). 



The energy of the complex amplitudes configuration for 
dipolar gas is given by: 
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where N is the total number of atoms and Ejnt {{cti}) is the 
quartic polynomial in the amplitudes {ai}. The sum is over 
all indices i, which obey inequality i^+iy < Kmax- The 
cut-off Kmax is the index of the last retained mode, defined 
as 2T:Kmax/L = k^ax- The interaction energies (dipolar and 
contact together) are given by 
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where n2D — N/L'^ is the two dimensional density. Here we 
introduced an effective contact interaction coupling constant 
9 ~ 9 + 9dd (cos 20 + cos^ 9) , where 9 is an angle between 
the direction of the dipoles' polarization and the Z axis. The 
second part of the dipolar interaction E^J^^ consists of the long 
range potential between dipoles: 
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with 



y^^(fe) = - 



2,1, kl cos 29 + kl cos^ 9 



e*^^'^/2Erfc ( ^ 
\/2 



(12) 

The derivation of formulas (10)-(12), with a short discussion 
of the potential (12) is presented in details in the Appendix A. 
The initial three dimensional problem is reduced to a two 
dimensional one. Hence the classical field 
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will refer from now on to a two dimensional function. 

From Eq. (10) it is clear that the dipole interaction in the 
geometry considered here partially acts like a contact interac- 
tion, thus it shifts the coupling constant g. Furthermore, this 
shift in the coupling constant may be tuned by the polarization 
direction. At the angle 9 « 0.96, at which cos 20+cos^ 6' = 0, 
this shift vanishes. 

At finite temperatures the probability distribution of a given 
configuration of the classical field according to the canonical 
ensemble is: 



Pi{o^^})^ 
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where Z{N, T) is the classical partition function for N atoms 
at temperature T. 

The average of any observable O is computed according to: 

{0)= /dao...dafc_P({aJ)0({aJ), (15) 
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FIG. 1. (Color online) Cuts through single exemplary copies of 
atomic density (three non-solid lines) obtained during generation of 
the canonical ensemble with the Metropolis algorithm, together with 
the average density over the ensemble (solid black line). 



where integrals go through all the possible complex ampli- 
tudes subjected to the constraint forced by the fixed number 
of particles: 
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= 1. 
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The best known Monte Carlo realization of this distribution is 
given by the Metropolis algorithm [23]. With it we explore 
the multidimensional phase space built from a set of {ai\ as 
described in detail in [17]. 

The algorithm realizes the idea of a statistical ensemble - 
it looks for suitable copies of the system (different configura- 
tions {a;}), all copies with energies distributed according to 
the Boltzmann distribution. As illustrated in Fig. 1 the single 
elements of the ensemble often differ significantly from the 
average of them. Typically, computations require millions of 
generated copies to reach reasonably smooth averages. 

Following [24], the identification of the condensate requires 
a diagonalization of the single particle density matrix: 
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where the mean value is taken with the probability distribution 
(14) and the eigenvector /3(0) corresponding to the dominant 
eigenvalue Ao (if such one exists) is the condensate. 

Let us show the outcome of the method for a special case 
of very low temperature. As the global phase has no physi- 
cal meaning, we can always eliminate ao choosing it as real 

and equal to ao = 



1-E 



fc#o Ffc"fcl 



The interaction en- 



ergy in (9) will remain a quartic polynomial in the amplitudes 
{aq^o}- If the amplitude of the zero momentum component 
is dominant, that is |ao| ~ 1 ^ lafc^^ol we can follow the 
derivation of the Bogoliubov approximation. It requires ne- 
glect of the third and the fourth order terms in the interaction 
energy. The remaining quadratic form is easily diagonalized 
yielding: 
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and the new complex amplitudes /3fe are defined as 
/3fc = Ukak - Vka*_f, 



(19) 
(20) 
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Thus the Bogoliubov spectrum, presented for instance in 
[25, 26], appears here as a special case. 



II. EXPERIMENTAL CONDITIONS 

Here we present the statistics of the quasi 2D dipolar gas 
for parameters close to the experimental ones, used in [1]. 
The experiment has been done in a deep optical lattice at tem- 
peratures around lOOnK. In each lattice site, the cloud was 
indeed quasi2D (checked with the help of numerically solved 
three dimensional Gross-Pitaevskii equation) with uJz/2tt = 
36573Hz. 

Our calculations are made for such parameters, that mimic 
the conditions in a single central lattice site. Thus we keep the 
aspect ratio between the Gaussian width in the Z direction and 
the size of the system in others equal to 100. To ensure a den- 
sity of atoms close to the experimental value (w 10^"' /cm'^) we 
consider here N = 500 atoms. We characterize the strength of 
the dipolar interaction with the "dipolar length" add ~ "l^fS 
[27], which for Chromium equals 16aB, where as is the Bohr 
radius. All results are presented in oscillatory units, i.e. length 
is given in l^ = \fhjrnc^z, where vicr is the Chromium 
mass. In a calculations for other atomic species, we use ^^^Er 



withafj = 70a f 



''dd ~ ' ""^s ^'^'1 ^^^Dy with a^J = 140aB, we keep the 
same density as for Chromium. Thus the geometry is changed 
in these cases: for atoms with mass m we choose the trapping 
frequency ui'^ = ^icr ^^ whereas the two dimensional box and 
the length unit h remain the same. 

For such parameters the fraction of atoms in the zero- 
momentum mode, shown in Fig. 2, is only slightly altered 
by interactions. We choose the scattering length IOqb, which 
is below "the dipolar length" add = lOas, to enhance the ef- 
fect of dipolar forces. We restrict the study to temperatures up 
to about lOOnK. Then the cutt-off in momentum space k^ax, 
limits also the theoretical spatial resolution. Here, the spatial 
resolution is > lOlz, which is not much smaller than the ob- 
served structure, which could slightly impact the outcome of 
the calculation. 

The dipolar interaction slows down a thermal depletion of 
the dominant mode. This is mostly because of dipole repul- 
sion - the experiment is in such a regime, that contribution of 
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FIG. 2. (Color online) Average fraction of condensed atoms. Pa- 
rameters: a = 10a_B, "the dipolar interaction length" corresponds 
to different atomic species ^^Cr, ^®*Er, ^^''Dy. The average density 
in the middle of the trap 10^* /cm''. Note, that oscillatory units are 
different for different species, as described in the text. The exact 
calculation for ideal gas follows the procedure from [22]. In inset: 
fraction of condensed Cr atoms at Ti — 0.015hujz/kB but for dif- 
ferent values of cut-off. 



the long range part of the dipolar potential (11) remains small. 
Thus, the main effect of the dipolar interaction is just a shift 
of the contact interaction coupling constant. It has already 
been proven [28], that contact interactions in the box reduce 
the evaporation to higher modes, as we also see in Fig. 2. 
The latter can be guessed also from the Bogoliubov spectrum 
(19): the energy gap between the and 1 momentum modes 
is increasing with both g and gdd- Note the opposite behav- 
ior for the gas trapped in a harmonic potential [12, 14]. The 
latter can be qualitatively understood with the help of the lo- 
cal density approximation, in which the system is considered 
as a set of boxes with varying chemical potential. If the lo- 
cal density nioc in a box decreases, while temperature is con- 
stant, then the fraction: average density of the atoms in the 
0-momentum mode divided by riioc, decreases also. Due to 
the repulsive forces the cloud trapped in a harmonic potential 
widens, and in consequence in the center of the trap (so in 
the region contributing mostly to the condensate) the density 
decreases. Thus in the dense part of the cloud atoms evap- 
orate to higher energy modes. The same reasoning leads to 
the conclusion, that in the tails of the cloud, where the local 
density would increase after the cloud broadening, the occu- 
pation of the 0-momentum mode increases. However, it turns 
out that globally, when repulsive forces increase, the conden- 
sate in trap evanesces. 

It is worth noticing, that the repulsive dipolar force has been 
already used to stabilize the system against attractive contact 
interactions [1] and even to suppress chemical reactions be- 
tween KRb molecules [29] and dipolar relaxation for spinor 
gases [30]. 

The discontinuities visible in Fig. 2 are due to the discret- 
ness of the finite-system spectrum. Precisely, there are regions 
where the temperature rises, and consequently cut-off energy, 
but the number of modes remains unchanged. At some value 
of temperature new four modes are accepted, which results 
in a small jump in statistical quantities. Our qualitative results 
do not depend strongly on this technical problem. Moreover in 
the rest of the paper we focus on such temperatures, at which 
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FIG. 3. (Color online) Cuts of the spatial first (the upper panel) and 
the third (the lower panel) order correlation functions at few chosen 
temperatures relevant to the experiment [1]. 



the energy of the last retained mode is equal exactly to the 
cut-off energy. The problem may be overcome by taking the 
weighted averages obtained for the "critical" cut-offs. 

In inset in Fig. 2 we investigated the sensitivity on the cut- 
off. We were changing the numerical factor in Eq. (8) keeping 
constant temperature Ti = 0.015fiujz/kB- Even for a cut-off 
varying up to 20% the fraction of condensed atoms is varying 
much less than 2%. One can see in the inset, that the statisti- 
cal eiTor is small - the main uncertainty comes from the dis- 
creetness of the spectrum, which cause the jump at the cut-off 
energy l.OS/ssT. 

The most interesting quantities are correlation functions. 
Among them, the crucial role for identifying condensation is 
played by the first order correlation function defined here as 



g\p)^{^*{p)yi>{0))/n2D. 



(21) 



In the 2D uniform system the true condensation does not ex- 
ist at r > 0. Indeed in Fig. 3a) one can see that even for 
very low temperatures g^{p) is decaying slightly below 1 for 
p 7^ 0. The width of this function defines the correlation 
length. It decreases with increasing temperature making the 
coherent region much smaller than the size of the system. In 
our regime, where both contact and dipolar interactions are 
weak (comparing to the contact interaction at the background 
scattering length), this reduction is mostly due to a shortening 
ofthe deBrogUe wavelength A(i = . '^ For temperatures 
used in Fig. 3 the A^ is varies between IIZ^ and 25lz- 

The possible density fluctuation are quantified with the help 
of the third order correlation function. Here we focus on its 




FIG. 4. (Color online) The first (the upper panel) and the third 
(the lower panel) order correlation functions at temperature T — 
OmOhu^/kB (^ 30nK for Cr ) for a = Was but for different 
values of add corresponding to different atomic species. 



versus scattering length is shown. For large scattering length 
the value of g^ does not depend strongly on the dipolar in- 
teractions, as it is insensitive to different polarization angles. 
This increase for smaller and smaller a is probably caused 
by approaching unstable regions. The latter is suggested by 
the rapid bunching of the least stable configuration 6 = tt/2. 
In such a geometry, when a reaches some critical value in the 
Bogoliubov spectrum imaginary frequencies appear, implying 
a phonon instability [33]. Then, in a collapse the gas would 
shrink almost to a delta function. Thus in the unstable region 
the system should be a mixture of infinitely many plane waves. 
This qualitative analysis indicates however that also close to 
instabilities, (but still in a stable regime) the occupation of 
fc 7^ modes should already increase. This is the regime 
studied throughout the paper - where one can see akeady dipo- 
lar effects, although the Bogoliubov spectrum does not reach 
roton yet. The macroscopically occupied higher modes give 
an increase of the g^ function. However, the dipolar interac- 
tions stabilize the system, due to their repulsive part. One can 
see already in Fig. 2 that the occupation of the 0-momentum 
mode increases with the increase of the dipole moment. Thus 
the role of excitations is getting smaller for stronger magnets, 
making the atomic cloud more coherent and less bunched, as 
shown in Fig 4. 
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FIG. 5. (Color online) The third order correlation function g^(0, 0). 
The role of dipolar interaction increases when approaching the insta- 
bility border. The temperature T — 0.02!iLOz/kB 
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We choose this quantity instead of g^ as it has a direct impact 
on the rate of three body losses [31], thus it is relevant for 
the condensate lifetime. Fig. 3b) shows, that surprisingly, in 
experiment [1] the gas was probably highly bunched, with a 
few times faster losses (at T > lOOnK) comparing with the 
pure condensate (T = OA). The bunching is simply due to 
approaching the thermal state, at which g"^ = 6 is expected. 

In Fig. 4 we compare the coiTelations of Chromium with 
the ones of Erbium and Dysprosium. As can be seen the cor- 
relation length is increasing with increasing dipolar forces and 
at the same time bunching is getting weaker [32]. These ef- 
fects may be understood with the help of Fig. 5 where g^{0, 0) 



III. CONCLUSIONS 



We have presented the classical field approximation for a 
weakly interacting dipolar gas. We show the statistical prop- 
erties relevant for present experiments with Chromium. The 
dipolar forces there are only a small perturbation around the 
ideal gas case. However, at the temperature at which the ex- 
periment [1] was done, the third order correlation function 
is highly increased compared to the zero temperature case. 
Due to this bunching, the three body losses should be a few 
times faster At the same temperature the correlation length 
is strongly reduced. In the case of atoms with stronger dipole 
moment. Erbium and Dysprosium, the reduction of the size of 
the coherent regions is getting smaller. 
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Appendix A: Derivation of the dipolar interaction energy 

Dipolar energy in the classical field approximation has a 
form: 

Ea.p = lJd^rJd^r'\^{r)\^\^ir')\^Vd,^{r^r'), 

(Al) 
where the three dimensional classical field 



K„ 



*3^(r) = VnMz) Y. 



gikip 



(A2) 



is a counterpart of the field operator and Vdip [r) is a dipo- 
lar interaction potential for atoms polarized in the Z direction 



ary conditions enforce discrete values of momenta (3). The 
integral j d^ r in (Al) is over all points from the domain 
[— ^, ■§■] X [— 1^, ■§] X R, where functions have Fourier trans- 
form of the form 

h{r) = ^^j ^Y.^i (fc.) 6*'='^+^'=^^ (A3) 



Assuming dipoles oriented within the XZ-plane the expan- 
sion coefficient of the dipolar potential is equal to 



Vi{k, 



J2^\ k^ + kj 



(A4) 



given in Eq. (2). Here i denotes a two-dimensional index where is an angle between polarization of all dipoles and the 
[ia;,iy] and the sum in (A2) is over all indices i which ful- 



fil the condition 
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i'i < Kmax- The periodic bound- 
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Combining Eqns. (A4), (A3) and (A2) in Eq. (Al) one gets 



Z direction and k^' = '^[ix, iy] is the x component of the 
vector k and fcj is the projection of fc onto XY - plane. 
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Thus the expression for the dipolar energy consists of a prod- for brevity we introduce an operator: 
uct of four separate integrals (related to Fourier transforms), 
afterwards integrated over k^ ■ There are two integrals which 
are Fourier transforms of i/)g : 
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i^l{z)e^'''^'Az = e'''"'"'\ 



(A6) 



Both integrals over coordinates in the box lead to the Kro- 
necker delta functions: 
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d'p J2 «:„a„e-*(— ±'^-)'' = L^J2 ^*m<^r.5ti^r 



Thus the formula for dipolar interaction energy (Al) re- 
duces to: 

^d-ip — 2L7 J y2¥ 

i,m,n,p,q'^m^p'^q^nOyi_ynOq_pVi (Kz ) 
— 2 J /9 — -^ mpn^ n—m \^z) i 

where in the last equation two delta functions reduce sums 
over i and q. With n2D — N jl? we denote a 2D density and 



Hence, the formula for dipolar interaction energy is reduced 
just to one integral (over fc^) and sum over three 2D indices 
(m, n, p). Substituting finally the explicit expression for the 
dipolar potential (A4) we have to deal with an integral 
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I^'Mx 
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sin^ Q + kl cos^ 6* 
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-1 
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The dipolar potential may be rearranged 

3^ -^ 1 = cos 26* + cos^ B 

kl cos 26 + kl cos^ 9 

to a form more convenient for integral over k^- This rear- 
rangement allows us to write the integral over momenta in the 
Z direction (A7) as a combination of two converging inte- 
grals: 
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FIG. 6. (Color online) a) Function V^^ defined as in the equation 
(AlO) for 6 = (polarization in Z direction) in momentum space, b) 
Function V^^ defined as in the equation (AlO) for 6 = (polariza- 
tion in Z direction) in position and c) its scaling for small and large 
P(c). 



where kn-m = '^yi'^x - n^)^ + {my - Uyf is just a 
length of the difference km, and fc„. Finally the dipolar en- 
ergy is a sum of two terms 

I Part of dipolar potential which mimics the contact in- 
teractions: 



E = 



9ddn2DN (cos 26 + cos^ 9) , 



2V2^l 



(A8) 



II Part of dipolar potential which depends on momenta in 
a non-trivial way: 
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3L kl cos 29 + kl cos^ 9 ,2 
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We illustrate the features of the potential in a particular case 
of 6' = 0. As shown in Fig. 6a) in this case the function (AlO) 
for fc ^- 00 converges to a non-zero value. In the position 
representation this asymptotic behavior leads to an additional 
delta function, precisely: 

T{v'^{k))=-\^^5{p) + V'''{pl (All) 

where the potential in a position representation is given by 



and 



*. 



-klll/2 



^[kr,-rr.l.r /2 



'27r k + fej; 



V ^fcn— T 



-Erfc 






F^^(p) = 27r / dfc Jo (M (l^'^ (fc) + y ^ ) fc' 



(A12) 
where Jq (x) is a zero order Bessel function. The shape of this 
potential is plotted in Fig. 6b). It diverges logarithmically at 
origin and tends to like ^ for large p, as shown in Fig. 6c). 
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